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We study distributions of eigenvalue curvatures for a block diagonal random matrix perturbed by 
a full random matrix. The most natural physical realization of this model is a quantum chaotic sys- 
tern with some inherent symmetry, such that its energy levels form two independent subsequences, 
subject to a generic perturbation which does not respect the symmetry. We describe analytically 
a crossover in the form of a curvature distribution with a tunable parameter namely the ratio of 
inter /intra subsystem coupling strengths. We find that the peak value of the curvature distribution 
is much more sensitive to the changes in this parameter than the power law tail behaviour. This 
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observation may help to clarify some qualitative features of the curvature distributions observed 
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experimentally in acoustic resonances of quartz blocks. 
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I. INTRODUCTION 



Studies of statistical properties of complex quantum systems (chaotic or disordered) show 
that their eigenvalue spectra exhibit patterns of universal fluctuations, whose structure 
mainly depends on the fundamental symmetries of the Hamiltonian Such a uni- 

versality opens an attractive possibility of modelling the fluctuations by comparing them 
with those observed in long sequences of eigenvalues of random matrices of appropriate 
symmetry ^, |^ . Namely, systems with no time reversal invariance are known to be ade- 
quately described by Gaussian Unitary Ensemble (GUE) of complex Hermitian matrices, 
and systems with time reversal invariance are described by Gaussian Orthogonal (GOE) or 
Gaussian Symplectic (GSE) ensembles of real symmetric or complex quaternion matrices, 
depending on the existence of strong spin-orbit coupling. 

More recently, the interest in studying spectral statistics of quantum chaotic systems was 
much revitalized by the understanding that their energy spectra display a universal change 
in their characteristics as a response to external perturbations of various kinds. The nature 
of the perturbation may vary considerably depending on the physical system, and usually 
involves application of external fields (magnetic or electric), change of boundary conditions 
or the shape of the system, rearrangement of positions of impurities in disordered medium, 
or variation of temperature, pressure or any other tunable physical characteristics. 

One of the frequently used measures of parametric sensitivities of complex quantum sys- 
tems is the distribution of level curvatures, which are defined as the second order derivative 
of eigenvalues with respect to a perturbation parameter. 

In Gaspard and co-authors developed expressions for the probability densities of level 
curvatures and found that the large curvatures must exhibit universal behaviour classified 
according to the underlying gross symmetries. Indeed, the curvatures become large in the 
vicinity of avoided crossings of energy levels as functions of a parameter and their distribu- 
tion can be simply related to that of small eigenvalue spacings. On the basis of extensive 
numerical investigations of both random matrices and several quantum chaotic systems, 
Zakrzewski and Delande conjectured an analytical expression for the full distribution of cur- 
vatures Later on Zakrzewski-Delande formulae were derived analytically by von Oppen 
[3] and by Fyodorov and Sommers [^] for the random matrix models of all three universality 
classes. 
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However, in many relevant experimental circumstances, physical systems have acciden- 
tally more underlying symmetries (frequently called "geometric") acting in addition to the 
presence or the absence of the time-reversal invariance. Such symmetries naturally induce 
classification of energy levels according to irreducible representations of the corresponding 
symmetry group, and the energy levels corresponding to different representations form sta- 
tistically independent subsequences. Only these subsequences may be then meaningfully 
compared with the universal random matrix patterns. In fact, a generic situation may be 
even more complicated, since geometric symmetries may not be exact, but approximate. 
This will clearly lead to spectra being a mixture of different subsequences with uncertain 
statistical consequences of mutual interference. 

Recent experimental studies on acoustic resonance spectra in quartz blocks 0,[l(| suggest 
that the system may fall into the latter category, and the deviations from standard theo- 
retical predictions of parametric correlations may have their origin in remnant geometric 
symmetries. This fact motivated several grou ps to investigate the effects of partial symme- 
try breaking on the level curvatures [ill ll^. However, closed form analytical expressions 
for curvature distributions for the case of partly broken symmetries are not available, to the 
best of our knowledge. 

In the present paper we consider a simpler, but related model. Rather than studying 
level curvatures in a system with partly broken symmetry, we address the case of two non- 
interacting subsystems subject to a perturbation which induces both coupling between the 
subsystems and variation of parameters within each of the two subsystems. This should be 
generically the case for a perturbation which does not respect the underlying symmetry. Of 
course, we do not claim that our simple model to be adequate for describing experimental 
situation in quartz blocks, but we rather hope that our results could help in indicating how 
various factors may affect the shape of level curvature distributions. 

Employing the random matrix calculations allows to derive exact expressions for the level 
curvature distribution as a function of relative weight of induced inter- and intra-sublattice 
osculating variations. The curvature distribution naturally interpolate between the simple 
Cauchy-Lorentz shape - characteristic for a pure symmetry-breaking perturbations which 
couple two subsystems without modifying them individually, - and the Zakrzewski-Delande 
formulas typical for perturbations respecting the underlying symmetry. 

We first treat a simpler case of complex Hermitian matrices (systems with broken time 
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reversal invariance) in detail, and then extend the derivations to the case of real symmetric 
matrices. Our analytical calculations are supported and corroborated by accurate numerical 
simulations of random matrix ensembles. 



II. GENERAL RELATIONS 

To study the parametric dependence of energy levels of a system with some underlying 
symmetry, we consider a random matrix model where the Hamiltonian 7i of the system 
linearly depends on a perturbation parameter e: 

H(e)=A + eB. (1) 

/ #1 

For the unperturbed Hamiltonian A we choose a block-diagonal matrix A = 

\ H 2 

with H\ 2 being N x N random matrices (complex Hermitian or real symmetric) taken 
either from GUE or, respectively, from GOE. These can be thought of as representing two 
non interacting chaotic subsystems. That means, we consider the matrices H p as random 
Gaussian, with entries H^j being independent and identically distributed variables with 
mean zero and variances < HijHy >= 2a 2 /N for the GUE case ((3 — 2), and < Hfj >= 
a 2 /2N for GOE case ((3 = 1). The joint probability density for H p is then written as 

N „ 



V(H p ) = CP N exp\-—TrH p \, (2) 



where C 1 ^ is the appropriate normalization constant. 

We introduce an interaction between blocks by considering a general coupling matrix 

I Ji v\ 

B = \ . We assume here that Ji 2 have the same symmetry properties as the 

\ V* J 2 J 

diagonal blocks H\ 2 and their probability densities can be written as 

P(J p )ocexp|-^TrJ p 2 |. (3) 

As for the off-diagonal blocks V, they represent general Gaussian random matrices (complex 
for (3 = 2 or real for (3 = 1) with no further symmetry constraints are imposed. The 
probability density for V is then chosen for both cases to be 

V(V) oc exp I ~^2 Tr ^y } ■ ( 4 ) 



Models of this kind were previously employed with satisfactory results in the analysis of 
data relative to symmetry breaking in nuclear physics Q|. 

Denote the eigenvalues and corresponding eigenvectors of H p by [A^, ] with p = 1,2. 
This implies that H p v^ = A^V^ where i — 1, . . . , N and v,- = 1. Our main goal is to 
find the distribution of level curvatures, defined as the second order derivative of eigenvalues 
of 7i with respect to the perturbation parameter e. Employing in the usual way the second 
order perturbation theory we can write the expression for the curvature corresponding, say, 
to an eigenvalue that for e = coincides with the (unperturbed) eigenvalue A^ as: 

r _ ( v Ii-A v ii)( v !i^i v fci) | (v!i^Vfc 2 )(vL^ tv *i) (r\ 

k^i A i A k k=l A i A k 

This expression shows that there are generically two contributions to the level curvatures. 
The first sum is essentially the level curvature induced by the block- diagonal part of the 
perturbation which does not lead to any mixing between levels of the two non-interacting 
subsystems. Taken alone, this term (which we will denote here as Cu) must, therefore, yield 
the Zakrzewski-Delande curvature distribution. In contrast, the second sum that will be 
denoted as C21 reflects the influence of the off-diagonal perturbation, which mixes the levels 
of two subsystems. 

Then the distribution of the total curvatures defined as V(C) = (5 (C — Cu — C21)) n and 
it can be conveniently written using the Fourier transform as 

V{C) = i- jH dxe- c {{e- c -) k (e- lxC %) M2 ■ (6) 

The first factor (^e~ lxCli ^ j J . is just the Fourier transform of the known Zakrzewski- 
Delande expression, hence the calculation of the curvature distribution reduces to evalu- 
ating the remaining factor (( e ~ lxC21 ) v ) ^ ■ Our next goal is to derive the corresponding 
expressions, first for (3 = 2 and then for /3 = 1. 



III. COMPLEX HERMITIAN MATRICES: (3 = 2 

We are going to evaluate the following ensemble average (see Eq.©) 

e A * - x k j . (7) 

/ V,H 2 



First, we perform the average over V. To simplify the notation we denote x/{X\ — X k ) = 
b k , v| 1 V r v fe2 = Wk and use the identity 

where the integration is taken over an auxiliary complex variable Z k . This allows to rewrite 
Eq-0 in the following form: 

Explicitly, employing the distribution function in Eq.(HJ), we need to calculate the integral 

dVdV^exp^-^Tr(VV r ) -iv^V ^^v fe2 j ^^fcvj 2 j ^ij • (10) 

It is convenient to introduce an N x N matrix G = (J2k ^k^k2) ® v Ji, so that, when used 
in the former identity relation we get 

dVdV^e oc e -^ Tr ( GtG ), (11) 



which is a generalization of Eq.fjHJ). We further use 

Tr(&d) = (v| lVil ) ^^^(vi 2 v, 2 ), (12) 

and recall that Vj! are eigenvectors of Hi and v k2 are those of H 2 . Using the orthogonality 
of the eigenvectors: 

1, k = g 
0, k ^ q 



and vjjVji = 1, gives the result for the average in Eq.tjHJ) as 



v*=i- 2 * / A \det 

For this we also performed the Gaussian integrations over Z k explicitly, and denoted x 

In what follows we denote A^ = A and Xj 1 + ^ = A& and proceed with calculations of 



the average (•••)# in Eq.(JT3J) by employing a technique suggested in |Tj|. In fact, for (3 = 2 
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he averages of the ratios of determinants are known in full generality for any value of iV 
13 , [l^ . Nevertheless, we outline the corresponding calculation in order to introduce the 

method and the convenient notation which will be used later on in this paper for the more 

complicated case [3 = 1. 

Using the standard "supersymmetrization" idea 16], we represent the denominator of 

the expression to be averaged as a Gaussian integral (we assume here x > for definiteness) 

det-\X b I N -H 2 ) = ± f dSd&e^ s -^"* s , (14) 

with a complex N dimensional vector S = (Si, . . . , Sn) t , where T stands for the vector 
transpose. For the determinant in the numerator, we use Gaussian integrals over anticom- 
muting (Grassmannian) N- component vectors \i X > which gives 

det(XI N -H 2 ) = L j d x d x t e *A X t x -! x t£ aX _ (15) 

Substituting the relations: x* H 2 x = ~Tr(H 2 x ® an d S^H 2 S = Tr(H 2 S (g) S^) in the 
integral, yields 

(■■■)h 2 = j d\j ^SgHAxt^sts]^-^^-^]^ (16) 

The ensemble average over GUE matrices H 2 can be easily performed by exploiting the 
identity 



e -iTr[HA]\ oce-^Tr^] (17) 

/ GUE 

and "decoupling" the quartic term in Grassmann variables with the help of simple Hubbard- 
Stratonovich transformation: 
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e m(x'x) 2 = f°° J2- e -4"Jfcx\ (18) 
J -oo V2~k 

After some straightforward manipulations we arrive at the following integral representa- 
tion for the required ensemble average: 



2 J J-oo \2ir 



1 ... Qcr a. S (g) S 1 " 
2 {lX ~Nj2 )lN ~2Np 



(19) 



Introducing the variable qp = iX — ^Jj^ an d shifting the contour of integration in such a 
way that, the integral over qp goes along the real axis, we can rewrite the above expression 
as 



(•■■>*« r^Le-^-^ 2 /d 2 Se^ sts ) 2 +^ sts det 



2 

q F I-?-S®& 



(20) 



where we shifted the contour for q F e (—00, 00) to be real. Further simplification can be 
made by noticing that the N x N matrix S ® SMs of rank unity, i.e. it has (N — 1) zero 
eigenvalues, and only one nonzero eigenvalue equal to (SS^). Then the determinant in the 
previous expression is equal to 

r 2 



det 



2 

q F I - ® S f 



q N F ~ l ( OF - ^SSt ) . (21) 



Finally, we introduce polar coordinates: S = rn with n^n = 1 and J d 2 S = r 2N ~ 1 dr dn, 
where J dn = produces a constant factor, which corresponds to the area of a 2N dimen- 
sional unit sphere. Further introducing p = r 2 and changing p — > Np/a 2 and then following 
with the obvious manipulations we get: 

/ det(AJjV - H) \ _ ( j Ne £ T X 2 f°° d( lF c -J^(gl-2i\ gF -^lnq F ) 



< ^ ) = U N e& n / —=^- 

\ det(X b I N - H) I ^ J -00 V2^qi 

x r dp(q F -p) c _ x 4 Pc _j^ (p 2_ 2iXp _ 4a 2 lnp ^ ^ 
Jo V2np 

where we reinstated A& = A + ix/N, x = 2a 2 x and stands for the accumulated constant 
factors. The latter can always be restored by noticing that when A;, = A the right hand side 
must yield unity identically. 

So far all the expressions were valid for finite-size matrices. When N — > 00 we expect the 
results, when appropriately scaled, to be universal, i.e. broadly insensitive to the details of 
the distribution of random matrices and applicable to quantum chaotic systems. In such a 
limit the integrals in Eq. (j22|) can be evaluated by the saddle point method. For A < V8cr 2 
( the so called bulk of the spectrum) the relevant saddle points are p s p = tA+ ^°" 2 ~ — and 
qS^p _ i\±Vsa 2 -> J^ ^ j g eag y ^ gee only the choice q s F p = tA ~ v/8 2 °" 2 " — yields the leading- 
order contribution, due to the presence of the factor (qp — p) in the integrand. Substituting 
this choice into the integrand in Eq. (j22j) and evaluating the Gaussian fluctuations around 
the saddle-point values, finally yields 



det(XI N -H) 
det[(X + t§)I N -H] 



x>0 



exp |-|[zA + V8a 2 - A 2 ] . (23) 



It is easy to repeat the calculation for x < and find that for any real value of x the result 
can be written as 

det(A/ iV - H) _\ 
HI / 



det[(\ + i§)I N -H] 



exp\-i\x^-2\x\np(\)a 2 }, (21) 



where p(A) = j^y/Sc 2 — A 2 is the mean eigenvalue density for GUE. 

The Fourier-transform of the above expression with respect to x immediately gives us the 
distribution V ff(C) of level curvatures induced by purely off-diagonal random coupling V 
between the two subsystems. In the large-size limit N — > oo we therefore have: 

-p (r s 1 Za 2 2 iT P (\) 

v off(C) = -- — 2 , (25) 



7T 



(C-A^r) + (2MA)o3)* 



<r 2 



/Ji 11 r. (i) ixi . (i)i ' 

A,- T7 A. 



which is nothing else but the Cauchy-Lorentz distribution with the mean value (C) a ff = A„ L > 
and characteristic widths T ff = (27rp(A)of). 

Turning our attention to the curvatures induced by the block-diagonal contributions 
J p (the term Cn in Eq.(jSJ) we can first perform the ensemble average over the Gaussian 
distribution of J\, Eq.fjHl). Employing similar methods as before, we easily find the result to 
be 

(Aj 1} -A ^ 

(1) ixi 

k^i y\ - — - * k 

where in this expression x± = 2o\x. This expression remains to be averaged over the joint 
probability density of N — 1 GUE eigenvalues \[ , . . . , A^i\, A^\, . . . , A^- which are different 
from the chosen eigenvalue A- whose curvature we address. The consideration which is 
exposed in 3, 0| shows that 

where H is a (N - 1) x (N - 1) GUE matrix. 

The averaging of the ratios of determinants in Eq.()27|) can be done, mutatis mutandis, 
by the same "supersymmetrization" procedure as above. The detailed exposition of the 
corresponding calculation can be found in the paper by Fyodorov and StrahovQ]. Here we 
briefly sketch the main steps. After representing each of four determinants by the Gaussian 
integrals (three over anticommuting and one over usual complex variables) one can easily 
perform the GUE average by exploiting the identity Eq. (jl7|) . Then the terms in the ex- 
ponent quartic with respect to anticommuting variables are "decoupled" by introducing an 
auxiliary integration over 3x3 Hermitian matrix Q, the procedure being a straightforward 
generalization of the Hubbard- Stratonovich transformation (|18j). All the subsequent manip- 
ulations are quite analogous to those exposed above, and for our case, instead of Eq.([22|) we 
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arrive at its analogue pertinent: 



/ det 3 (A H) \ ^ f dQ {detQF f~i exp(--^Tr(Q F -i\lf 
\det(X b -H)/~ J ( 4a 2 



x J~dpp*- 1 e-'$*exp I -£- 2 (p 2 - 2z\p) I -p) (<#> - p) (q® - p) , (28) 



where q F ' are real eigenvalues of the Hermitian 3x3 matrix Q p and N stands for N — l. In 
fact, since Tr(Qp — iXI) 2 and det Qp depend only on the eigenvalues qp' 2 ' , it is convenient 
to use these eigenvalues and the corresponding eigenvectors as integration variables. In these 
coordinates the integration measure is given by 

dQ F oc d^dqfdqfdqf \{ (qf ] - qf ] ) , 

l<fei<fe 2 <3 

where dfi[U] is the invariant measure on the manifold of unitary 3x3 matrices, representing 
the eigenvectors of Qp and the last factor is the Jacobian of the transformation, known as 
the Vandermonde determinant. 

Again, we are interested in the limit N 3> 1, so we neglect the difference between N and 
iV — 1 and omit the tilde henceforth. The set of the saddle points of the integrand with 
respect to each of the variables p > and g^' 2,3) is: p s - p = iX +^f~* an d q s / = ^f^ . 
These saddle points are the same as what we found earlier. However, the presence of both the 
Vandermonde factors and that of the product Hfc=i (<?f ^ — p) m& ke us select the following 
saddle points: 

(1) iX + V8a 2 - X 2 (2) (3) iX - V8a 2 - A 2 
Qf = 2 ' Qf = Qf = 2 

(as well as its cyclic permutations) as these give the leading-order contribution. In fact, 
the integrand vanishes at these saddle-point values and care should be taken to expand the 
integrand further when calculating the contribution from the Gaussian fluctuations around 
the saddle-points (see |3] for a general procedure). The final result is given by 



( e ~ lXCvi ) AA = I 1 + 2 MAK 2 M] exp {-*Ax|l - 2|xMA)a 2 } . 



(29) 



Taking the Fourier-transform, we, as expected, arrive at the Zakrzewski-Delande formula 
for (3 = 2: 

V diag {C) = 3- (30) 

*[(C-(C} d ) 2 + T 2 } 2 
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where the mean value (C} d = X-p? and the characteristic widths T d = (2tt p(X)af). 

Now we know all the factors in Eq.© and can find the curvature distribution accounting 
for both diagonal and off-diagonal perturbations of the two decoupled subsystems: 

r°° dr f A 1 

V(C) = J ^ -e» c [1 + T d \x\] exp [-^(^ + ~ H ^ff + ■ ( 31 ) 

Performing the integration explicitly, we arrive at our final formula for complex Hermitian 
case: 

V(C) = - \ (32) 



2T d (T off + T 



*2 



d 



(C-l((C) off + (C) d )) 2 + (T off + T d ) 

IV. REAL SYMMETRIC MATRICES: [5 = 1 

We again need to evaluate the ensemble average as in Eq.(|7|). but this time for the real- 
valued perturbation V and real- valued eigenvectors Vj, so that the quantity vf 1 Vvk2 = Wk 
is a real variable. As before we denote x/(X[ — ) = bk and use the integration over an 
auxiliary real variable x^' 

-ib k wl _ . 1 f -±-x\-ix k w k 



V h J 2tx 
combined with the fact that (cf. Eqs. (jTTjl ) 

J dVexp ^-^L T r(VV T ) - iTrvfy k (v fc2 ® vg) J oc exp ^x 2 \ (34) 

because of the orthogonality of eigenvectors. Consequently, we easily perform the Gaussian 
integral over Xk and obtain 

/ det 1 / 2 (\fh N - H 2 ) \ 



det^ [( A ?> + §) ^ - # 2 



H 2 



as before we denoted x = 2o\x. 



After denoting A,- = A and \\ + S = A& for a less cumbersome expression, we then 
proceed with calculations of the average (•••)_#■ in Eq.(|35p. To be able to employ the 
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previous technique for (3 = 2 case we first rewrite: 

(let 1 2 (xi N - H 2 ) det ( XI y - H 



det 1/2 [ X h I N - Ho) det 1/2 ( X b I N - H 2 ) det 1/2 ( XI N - H- 



(36) 



Assuming, for definiteness, x < 0, and also assuming that A has an infinitesimal negative 
imaginary part we can represent the two factors in the denominator as Gaussian integrals 
over real N- component vectors Xi 2 : 



and 



det- 1 / 2 (A/ JV - H 2 ) oc J d Xl e-i Ax ^ Xl+ ^ x ^^ Xl 



det- 1/2 (X b I N -H 2 ) oc / c / X2 e-f A " x 2X2+fx 2 H 2 x 2 



(37) 



(38) 



where T stands for vector transpose. As for the determinant in the numerator, we can 
use the same Gaussian integral Eq. ffTo^) over anticommuting (Grassmannian) N- component 
vectors Xi X^ ■ Substituting these integral representations to Eq. (f3T)|) and performing the 
ensemble averaging over GOE matrix H2, with the help of identity: 



e ±iTr[M]\ K e -^Tr(A^+A)\ 

I GOE 



one can satisfy oneself that the resulting expression takes the form: 

dxidx 2 dx d X te-K^i+^i^-Axtx) 



(39) 



(40) 



x exp 



a 



IN' 



Tr 



Q 2 



2 2 

+ ^ (X j x) 2 + I ^X t (xi ® xf + x 2 ® x^) X 



T T 
X^Xi x{x 2 



T T 
X2 X]^ X2 X 2 



Now 



In this expression we introduced a positive definite matrix Q 

we again use the "decoupling" of the quartic term in Grassmann variables (the simple 
Hubbard-Stratonovich transformation Eq. (jl8|0 and then perform the Gaussian integration 
over anticommuting variables explicitly. The latter yields the determinant factor 

qa 



det 



iX- 



o- 



X ) In + 2N ( Xl ® X ^ + X2 ® X ^ 



(41) 



This factor can be brought to a simpler form 

N-2 



iX- 



qa 



det 



iX- 



qa 



a 

2N 



T T 

xfxi x|x 2 

T T 
XgXi X 2 X 2 



(42) 
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by noticing that (xi <8> + x 2 <8> x^) = XX T , where X = (xi,x 2 ) is iV x 2 rectangular 
matrix, and using the identity: det (I N - XX T ) = det (I 2 - X T X) then recognizing that 
Q introduced by us above is just 2x2 matrix X T X. We see that the resulting expression 
depends on the vectors X1.2 only via the matrix Q. In recent papers it was shown that 
the integration over x 12 under these conditions can be replaced by that over Q, with an 
extra factor detQ^ N ~ 3 ^ 2 arising in the integration measure. After some straightforward 
manipulations we arrive at the following integral representation for the required ensemble 
average: 

POO / \ N — 2 n 

(• dqe-f (iX--%) I dQdetQ^' 2 



x exp 



a 



8N' 



■Tr 



Q 1 



Tr 

2 



Q 



X 
Xh 



'N 
det 



iX- 



Q>0 

go- 



er 



V , h + 2N Q 



(43) 



Introducing the variables qp = —iX + and Q b = and shifting the contour of 

integration in such a way that the integral over qp goes along the real axis, we can rewrite 
the above expression as 



dq F qp 2 e 



r{q F +i\) 2 



x exp 



6>0 

A' 



2(T 2 



-Tr 



dQ b det Q b {N ~ 3)/2 det \-q F i 2 + Q 
Q 



iN 

—Tr 




(44) 



At the next step we introduce appropriate polar coordinates in the space of matrices Q b > 0: 



Q b = O t \ Pl ° | O , dQ b oc \ Pl - P2 \d Pl dp 2 dO 

P 2 



(45) 



where p± >2 > and O are 2x2 real orthogonal matrices: O O = I 2 , with dO being the 

(cos0 sin0 \ 
and 
— sin cos <fi J 

dO = dej) /(2ir). Substituting these expressions into Eq.(JSJ) and after obvious manipulations 

we get 

det 1/2 (XI N - H) \ = n X 2 f°° dqp N (q 2 F+2iXqF _ a 2 lnqF) 
det^(X b i N -H)/ &2 ' ° 



00 Qf 



x /„ **/„ dp2 M>» 



(qF-Pl)(qF-P2)Tx(Pl,P2)e 2 = 



r 0C(pi)+£(P2)) 



(46) 
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where we reinstated A& = A + ix/N, x = 2<j\x and 



C(p) = p 2 + 2i\p - a 2 hip , X x ( Pl ,p 2 ) = J ^ e 4 (pi+P2) - (pi - p2)cos20] 



and, as before, Cn stands for the accumulated constant factors. So far all expressions 
were valid for finite-size matrices. When N — > oo the integrals in Eq. (j46p are evaluated 



by the saddle point method. For the bulk of the GOE spectrum A < V2cr 2 , the relevant 
saddle points are p$ = an d q 8 / = ~W2^-a^ Again only the choice q 8 / = 

iA+V2<r 2 - A £ yi e i(j s leading-order contribution, due to presence of the factors (qp — Pi.,2) 
in the integrand. Substituting this choice into the integrand in Eq. (|4fj)l and evaluating the 
Gaussian fluctuations around the saddle-point values finally yields 

det 1/2 (A/ JV - J ff) 



det l / 2 [(\ + i§)I N -H] 



x<0 
N— >oo 



exp |x[-iA + V2a 2 - A 2 ] ■ (47) 



It is easy to repeat the calculation for x > and find that for any real value of x the result 
can be written as 



det 1/2 (XI N -H) 
det 1/2 [(A + if 



exp |— iXx— | — |a;|7rp(A)cr 2 | , (48) 



where p(A) = ^v^er 2 — A 2 is the mean eigenvalue density for GOE. 

We see that in the large-size limit N — > oo the expression for the level curvature distri- 
bution induced by purely off-diagonal random coupling V between the two subsystems is 
essentially the same Cauchy-Lorentz distribution for both (3 = 2 and [3 = 1 cases, up to 
re-scaling of the widths and the mean value with a simple factor 2: 

r«m = ~ ^ • m 



7T 



(C-\$) + (7rp(A)(T 2 ) 2 



We give a more detailed discussion of this issue in the next section. 

The distribution of curvatures induced by the block-diagonal contributions J p (the term 
Cu in 03) for real symmetric matrices is quite different from that of complex Hermitian ones. 
Performing the ensemble average over the Gaussian distribution of Ji, Eq.Q and employing 
the same methods we find the result to be 

- a (1 V/2 
[aJ 1} + ^ - a1 1} 
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with xi = 2a\x. The averaging over the joint probability density of (N — l) GOE eigenvalues 
X[ , . . . , A^i\, A^, . . ., Aw' which are different from the eigenvalue whose curvature we 
address, shows that p, |8[ 

where H is a (iV - 1) x (N - 1) GOE matrix. 

The averaging of the ratios of determinants in Eq. (j51j) can be done, mutatis mutandis, by 
the same technique as above. However, the presence of the absolute value of the determinant 
makes accurate calculation to be quite lengthy, and it will be presented elsewhere, but the 
result is compact and it is given by |3? Q] 

(e- ixCli )j u6l = vp^WMe-^K, (MMAK 2 ) , (52) 

with Ki(z) being the MacDonald function of the order one. Such an expression yields, after 
the Fourier-transform, the Zakrzewski-Delande formula for (3=1: 

1 r 2 

VdiagiC) = —jfi, (53) 

2 

where the mean value (C)d = — A^ and the characteristic widths Yd = np{X)a\. 

The curvature distribution, accounting for both the diagonal and the off-diagonal pertur- 
bations of two decoupled subsystems, can be found as the convolution of two distributions: 

/oo 
dC 1 P diag (C 1 )P ff(C-C 1 ) (54) 
-oo 

and in this way we arrive at the final formula for the real symmetric case. We present it 
below for the central point of the spectrum A = 0: 

P(C) = — / dCl -^-rj, . (55) 

2W-0C [c 1 2 + r 2 ] 3/2 (c-c 1 ) 2 + r 2 // 

V. NUMERICAL RESULTS AND DISCUSSIONS 

In the present section we compare the derived analytical form of the curvature distri- 
bution with the results of direct numerical simulations of the ensemble. For our numerical 
investigations we used a normal random distribution that was adopted from FORTRAN 
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Numerical Recipes |l7j and to find the eigenvalues we superseded some subroutines from 
LAPACK To avoid the necessity of unfolding the spectra we took into account only 
levels around the central part of the spectrum. Namely, for a 100 x 100 matrix, ten mid- 
dle eigenvalues (20 or more for larger matrices) were considered at each time step, and a 
curvature value for each eigenvalue was calculated by a second difference equation 



, , Ai(-2e) + 16Ai(-e) - 30^(0) + 16Ai(e) - A<(2e) 

\ { £ > \e=0 - " \ ' 



The choice of five points instead of the usual three jlj] was made to ensure the stability of the 
results, especially for the GOE case of our system. The empirical choice of e = 0.001 was an 
outcome of a number of trials; the values it takes may be system specific. We finally remark 
that using larger matrices e.g. 400 x 400 did not improve the quality of plots considerably. 

The normalized results of the simulations are presented in Figs. (|l|2jl . To compare them 
with the analytical predictions, for the GUE-like case (3 = 2 we consider Eq.(|32|) at A = 
and a = l. It is also convenient to use the dimensionless curvatures k obtained from C by 
reseat the latter with the variance of "level velocity as (of. fl)= 

« = C . A - „. (57) 




where A = airy/2 /N is the mean level spacing of a single subsystem at A = 0. After some 

Simp ,e calcalatioa, the , 8 t o.er pe—oa theo, gl ve S - ^ aa d 

is independent of inter-subsystem coupling strength cr 2 . As a result of the rescaling, the 
curvature distribution acquires the form 

?(«) = -(,, f m+ (1 + r)2 ,|. (58) 
w ir \ [<c 2 + (1 + r) ! ] [«2 + (1 + r ff \ 

controlled by the only parameter r = o\ja\ i.e the ratio of the inter-subsystem to the 
intra-subsystem coupling strengths. Superimposing the plots of this expression over the 
appropriately normalized numerical data shows good agreement for all corresponding values 
of the parameter r, despite the noise in the large curvature tails. For curvatures exceeding 
the typical value k ^> r + 1 the distribution shows a power law tail, the GUE-like behaviour 
k~ 4 being replaced by the Cauchy-Lorentz one k~ 2 with the relative growth of the ratio r. 
For any r > the most distant tail is always of the Cauchy-Lorentz type, but intermediate 
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GUE-like behaviour is clearly seen for r C 1 when intra-subsystem coupling appreciably 
exceeds the inter-subsystem one. The crossover curvature value between the two regimes of 
decay is approximately described by the expression 



(59) 



which can be obtained by equating the large-curvature tails originating from the two com- 
peting terms in the expression Eq.(|58p. For curvatures in the interval 1 <C «; C k ct ~ \/2/r, 



the behaviour is GUE-like, changing to a slower Cauchy-Lorentz decay at k> n cr ~ a/2/t. 
It is also easy to verify that the maximal value of the distribution P max = V (0) always 
decreases with the increasing ratio 



7T 



1 1 

+ 



(60) 



r + 1 ' (r + l) 2 _ 

For the case of GOE, the similar rescaling of curvatures at A = and a = \[2 leads from 
Eq. (jB3j) to the expression: 

1 f°° 1 r 

= 27 J ^ dX (x? + l) 3 / 2 r 2 + (x - Kf (61) 

We plot this distribution superimposed over the numerical data for various values of r. 
Again, they agree rather well with the numerics and a crossover behaviour from GOE-like 
tail to Cauchy-Lorentz one can be seen clearly for small r, e.g. r = 0.05. 

In fact, one may notice from our plots that decrease in maximal value of the distribution 
starts to be noticeable at much smaller values of r than the modification of tail behaviour 
at not very large values of k. This fact qualitatively corroborates with the experimental 
observations in quartz blocks where noticeable deviations were detected in the center of 
the distribution, whereas the tails agreed well. Although, our oversimplified model clearly 
can not be considered as adequate for describing the actual experimental situation, we 
nevertheless mention that the choice of r ~ 0.2 allows matching the drop of the peak value 
with the experimentally observed deviation and produces an overall good agreement with 
the experimental curve. 

As we already noted in the text of the paper, and as clearly seen from the numerical Log- 
Log plots, the limiting large —N curvature distribution due to purely off-diagonal (inter- 
subsystem) perturbations turned out to have the same Cauchy-Lorentz form irrespective 
of the underlying symmetry, for both (3 = 2 and (3=1. Below we give an alternative, 
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heuristic derivation of this fact, which sheds some light on the origin of such a behaviour. 
The starting point for our analysis is expression Eq.([13p for (3 = 2 or Eq. (j35|) for (3=1. 



Denoting A 



(i) 



A, as in the text above, we rewrite those formulas as: 
(A - A a , fc ) 



exp 



P 



N 



E 1 ^ 



k=l 



1 



N 

n 

\k=l 
IX 



(A+f) 

1 



N (A + f ) - X 2 , k 




IX 



a N — 
(3 x - ix 

~9 2-^ 



N (A + |) - X 2 , k 
1 



^(A + |)-A 2 , fe 




At the last step we made a plausible assumption that the expression above in the limit of 
large N can be approximated by expanding the logarithms in the exponential to the first 
non- vanishing term. Now we introduce an exact eigenvalue density for H 2 as: 

1 N 



k=l 



where 5(x) stands for the Dirac delta function. In terms of this density the ensemble average 
(|62[) can be rewritten as: 



exp 



-x-x 




(64) 



Now, we use the well-known fact that the exact eigenvalue density for random matrices is 
self- averaging, which means in the limit iV — > oo converges to a non-random smooth function 
-the mean eigenvalue density. The latter function is just given by the Wigner semicircular 
law jO sc (/i) = — V \J A*s C — A^ 2 for < /x sc , where fM sc = V8cr 2 for GUE and ji sc = \ / 2o~ 2 for 
GOE. All these facts suggest that in the limit of large N the ensemble average in (jfi4j) can 
be suppressed in favour of replacing the exact density with its semicircular form. Moreover, 
since p sc (/i) is a smooth function, in the limit of iV — > oo we can use the Sohotsky formula: 



lim 

N—fOc 



^ c (A + f ) - /i 



V 



- isgn[x\np sc (X), 



(65) 



where the first term is understood as a principal value integral, and sgn stands for the sign 
function of the argument. In fact, with some effort the integral can be evaluated explicitly: 

72 



V 



Msc ( A - N Use 



(66) 
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Collecting all terms we see that: 
det(XI N -H) 



det[(X + t§)I N 



H] 



13/2 \ 



-Pi"- 



.~2A 

ix— + x 7rp sc (AJ 

. /4 



(67) 



which coincides with the earlier derived expressions in Eq.([24j) and in Eq. (j48|) . 

It is natural to expect that such a derivation can be made mathematically rigorous. 
However, despite its simplicity and conceptual clarity, such a method can not be straightfor- 
wardly applied to evaluation of the more complicated averages such as those in Eqs.(J27J) and 
(I51J1 . Indeed, the application of the outlined procedure to Eq. (j27|) amounts to approximating 
the extra determinant factor in the large-N limit as: 

det 2 (\i N -H 2 ^ ^exp^N-^V dp yj ^ - p 2 log (A - p) \ 



e mm 



which is indeed a correct expression up to the leading order in N in the exponential. It 
serves to cancel the extra factor e~^ T in front of the ensemble average in Eq.([27|). However, 
it is easy to understand that to arrive to the correct expression Eq. (|29jl one needs to take 
into account subleading terms -those of the order of unity in the exponential. This goal goes 
beyond the simple use of self- averaging, and requires a much more detailed treatment. It is 
not clear at the moment how to implement such a treatment in the present heuristic scheme. 
That is why the "supersymmetrization" method, which yields fully controllable results in 
all cases should be in general preferred. 

In conclusion, we have derived exact expressions for the distribution of level curvatures in 
a model describing a mixing of two independent spectra by a generic perturbation. Although 
the model is too simple to describe actual experimental situation in systems with partially 
broken symmetries, some features of the behaviour of our curvature distribution may play 
a role of useful analogy helping to understand the deviations in experimentally measured 
level curvature distribution of the acoustic resonances of quartz blocks Indeed, the 
maximum value of the latter distribution was found to be considerably lower than predicted 
for pure GOE case, whereas the tail shows good C~ 3 decay. This agrees qualitatively with 
our observation that the peak value of the curvature distribution might be more sensitive 
to remnant symmetries than the power law tail behaviour. 

In fact, an ideal experimental realisation of our model may be the system of two super- 
conduction microwave billiards coupled by an antenna in a variable way 3]. Although, in 
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real experiments of this type the coupling was changed in large discrete increments, it is in 
principle possible to change it in a much more controllable way, and to study level dynamics 
induced by such a coupling. We hope that our results may stimulate experiments of this 
sort. 
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FIG. 1: Normalized GUE curvature distributions for a few selected values of r = o\ja\ (inter- 
coupling to intra-coupling strengths); bottom plot (Log-Log scale) shows the tail behaviour of the 
same distribution. 
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FIG. 2: Normalized GOE curvature distributions for a few selected values of r = (y\ja\. The solid 
lines are analytical predictions superimposed over numerical data plots. 
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